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conditions  are  first  explored  by  means  of  a  semidiscretization. 

(The  differential  operators  are  discretized  with  respect  to  the 
direction  normal  to  the  streamlines,  but  not  with  respect  to  the 
streamline  direction.)  The  resulting  system  of  ordinary  differen¬ 
tial  equations  has  singular  points,  whose  position  is  related 
to  the  transition  from  one  regime  to  the  other.  These  singularitie 
are  compatible  with  weak  solutions  of  the  problems.  (Finite 
difference  and  finite  element  solutions  are  realizations  of  the 
concept  of  weak  equality.)  But  from  the  physical  point  of  view 
these  singularities  are  not  admissible,  at  least  not  in  the 
transition  from  a  subsonic  to  a  supersonic  speed.  One  thus  obtains 
the  requirement  that  at  the  sonic  line  the  partial  differential 
equation  be  satisfied  in  the  strong  sense.  s  ^gl ves._addi t iona  1 

conditions,  which  are  first  expressed  in  the  semidiscretized  versio 
of  the  problem  and  later  introduced  directly  into  the  fully 
discretized  finite  element  version.  The  condition  so  obtained 
|  is  the  analogue  to  Murman ' s  sonic  operator.  The  singularities 
mentioned  above  cannot  be  avoided  in  the  return  from  supersonic 
to  subsonic  speeds.  Here  the  requirement  of  a  weak  solution  is 
appropriate.  Admitting  such  singularities  one  introduces  an 
|  arbitrariness  into  the  problem  which  is  necessary  to  satisfy  the 
boundary  conditions  at  the  subsonic  downstream  boundary  of  the 
*  flow  field.  The  admission  of  suitable  weak  solutions  is  equivalent 
to  shock  conditions.  In  the  present  context  it  is  required  that 
the  weights  be  constant  across  the  shocks.  Some  adjustment  of 
the  weight  areas  is  needed  in  order  to  match  the  number  of 
unknowns  with  the  number  of  equations.  The  shock  conditions  so 
obtained  are  correct  only  in  the  average  over  several  elements. 

This  might  lead  to  a  loss  of  accuracy  unless  one  chooses  the 
elements  sufficiently  small.  A  method  is  sketched  which  removes 
this  restriction  by  taking  the  shock  position  within  the  elements 
into  account.  A  formulation  of  this  kind  will  most  probably  be 
applied  in  combination  with  a  Newton-Raphson  procedure,  in  which 
one  determines  simultaneously  by  direct  elimination  corrections 
for  the  entire  flow  field  including  the  position  of  the  shock 
from  a  system  of  linear  equations. 
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SECTION  I 
INTRODUCTION 

The  computation  of  transonic  flow  fields  by  means  of 
finite  difference  methods  has  reached  a  high  degree  of 
sophistication.  One  of  the  most  successful  approaches 
(that  of  Jameson,  Ref.  1)  alternates  between  a  Poisson  solver 
and  an  iteration  process  patterned  after  that  of  Murman  and 
Cole  (Ref.  2,  3).  The  Poisson  solver  can  hardly  be  improved 
by  the  finite  element  concept,  but  there  is  a  chance  that  the 
task  performed  by  the  Murman  Cole  iteration  can  be  carried 
out  more  efficiently.  The  present  report  does  not  claim  to 
have  achieved  this  goal.  It  asks  a  preliminary  question: 
which  special  measures  are  needed  in  order  to  obtain  such  a 
finite  element  formulation.  The  approach  taken  is  primarily 
theoretical.  The  procedures  which  emerge  are  more  complicated 
than  those  of  Murman  and  Cole.  Whether  simplifications  are 
possible  remains  to  be  seen.  The  author  considers  it 
preferable  first  to  study  what  should  be  done  in  principle 
and  to  make  simplifications  and  compromises  only  after  the 
basic  procedure  has  been  established,  although  a  less 
systematic  approach  may  sometimes  be  equally  successful. 

There  is  no  doubt  that  a  finite  element  approach  is 
applicable  in  subsonic  problems,  for  the  computation  can  be 
based  on  an  extremum  principle.  Such  a  secure  basis  does 
not  exist  in  the  supersonic  region,  but  the  study  of  examples 
shows  that  stable  methods  using  a  finite  element  concept 
exist  (Ref.  4).  Two  approaches  are  possible:  one  can  either 
treat  the  subsonic  and  the  supersonic  regions  separately  and 
then  impose  matching  conditions,  or  one  can  try  to  develop 
a  procedure  in  which  the  subsonic  and  the  supersonic  regions 
are  treated  simultaneously,  although  not  exactly  in  the  same 
manner.  Even  in  such  an  approach  the  transition  from  the 


subsonic  to  the  supersonic  region  and  back  requires  special 
measures.  In  finite  difference  form  this  is  the  approach 
of  Murman  and  Cole.  In  the  present  report  the  analogous 
finite  element  version  will  be  studied.  Of  particular 
importance  are  the  conditions  that  must  be  imposed  at  the 
shock  and  at  the  sonic  line.  The  sonic  operator  and  the 
shock  operator  of  the  Murman-Cole  procedure  can  be  inter¬ 
preted  as  the  finite  difference  realization  of  the  ideas 
that  will  be  developed. 
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SECTION  II 


CONDITION  AT  THE  SONIC  LINE  OBTAINED 
BY  A  SEMI-DISCRETIZATION 

In  Reference  5  the  author  has  shown  that  a  special 
condition  must  be  introduced  at  the  sonic  line.  We  ask 
here  how  this  condition  expresses  itself  in  a  finite  element 
setting.  Reference  5  is  motivated  by  the  observation  that 
the  boundary  conditions  for  the  flow  in  a  channel  are  quite 
different  in  a  purely  supersonic  and  in  a  purely  subsonic 
flow.  In  a  supersonic  flow  one  prescribes  the  full  velocity 
vector  in  the  entrance  cross  section.  In  a  subsonic  flow 
one  has  conditions  in  the  entrance  and  in  the  exit  cross 
section.  If  one  has  a  transition  from  a  subsonic  to  super¬ 
sonic  flow,  one  has  subsonic  boundary  conditions  in  the 
entrance  and  no  boundary  conditions  in  the  exit  cross  section. 

In  some  manner,  one  must  have  extra  conditions  within  the 
flow  field  which  are  substitutes  for  the  conditions  which  would 
be  given  in  the  entrance  cross  section  if  the  flow  were  entirely 
supersonic.  The  nature  of  these  conditions  is  recognized  if 
one  treats  the  problem  by  means  of  a  semidiscretization;  the 
differential  operator  for  the  direction  normal  to  the  stream¬ 
lines  is  discretized,  but  the  derivatives  in  the  streamwise 
direction  are  retained.  One  obtains  a  system  of  ordinary 
differential  equations.  To  be  specific,  let  us  assume  that 
the  flow  field  is  described  by  the  values  of  the  potential 
at  a  discrete  set  of  streamlines.  (Mostly  we  think  of  the 
simplified  form  of  the  potential  equation  then  the  stream¬ 
lines  are  lines  y  =  const.  For  the  full  potential  equation, 
one  will  introduce  a  linearization  and  the  streamlines  are 
those  of  the  preceding  iteration  step.)  One  then  obtains  a 
system  of  ordinary  differential  equations  for  the  potential 
along  these  selected  streamlines.  One  finds  that  singular 
points  arise  where  these  lines  intersect  the  sonic  line. 


At  these  singular  points  all  but  one  of  the  solutions  of  the 
homogeneous  system  are  regular.  The  singularity  of  this 
exceptional  solution  is  compatible  with  a  weak  solution  of 
the  system  of  equations,  but  it  is  not  admissible  because 
then  the  streamwise  component  of  the  velocity  gradient  becomes 
infinite.  The  missing  condition  for  the  system  of  differential 
equations  is  obtained  by  the  requirement  that  such  behavior 
of  the  solution  is  excluded.  In  mathematical  terms  one  would 
say,  that  at  the  sonic  line  the  differential  equation  must  be 
satisfied  in  the  strong  sense.  For  the  nonlinearized 
differential  equation  this  amounts  to  the  requirement  that  the 
second  derivative  normal  to  the  streamline  of  the  potential 
be  zero.  This  is  the  condition  imposed  by  Murman.  For  the 
linearized  equation  which  one  must  solve  when  one  applies  a 
Newton  Raphson  procedure,  the  condition  is  slightly  more 
complicated.  It  is  best  derived  directly  from  the  linearized 
differential  equation.  The  following  analysis  differs 
somewhat  from  that  of  Reference  5  because  it  refers  to  a 
partial  differential  equation  which  resembles  more  closely 
the  one  which  arises  by  a  linearization  of  the  original  non¬ 
linear  problem. 

We  shall  deal  here  with  either  bilinear  or  bicubic 
shape  functions  in  a  rectangular  grid  system.  For  functions 
of  this  character  one  can  carry  out  the  discretization  for 
the  direction  normal  to  the  streamlines  and  for  the  streamline 
direction  separately.  Then  one  obtains  as  an  intermediate 
form  the  system  of  ordinary  differential  equations  mentioned 
above.  In  practice  this  intermediate  step  will  be  omitted, 
but  it  is  worthwhile  to  discuss  it  because  it  suggests  the 
form  of  the  final  numerical  formulation.  The  description 
of  the  practical  procedure  which  will  be  given  in  Section  III 
is  selfcontained.  A  reader  primarily  interested  in  this 
aspect  may  pass  over  the  present  discussions,  whose  main 
purpose  is  to  provide  a  motivation. 
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We  start  from  the  simplified  equation  for  transonic 


flow 

-(Y+l)Vxx  +  *yy  =  0 

Assume  that  one  has  an  approximation 

<i>(x,y)  =  F(x,y) 

and  that  one  tries  to  obtain  corrections  by  means  of  a  Newton 
Raphson  procedure.  Then  one  sets 

<J>(x,y)  =  F(x,y)  +  <f>(x,y) 

and  obtains  by  a  linearization 

-(y+1)  [F  d>  +  F  <b„]  <t>  -(y+l)F  F  +  F  =0 

w  '  L  xyxx  xxyxJ  yy y  w  '  x  xx  yy 

Set 

G(x,y)  =  -(Y+l)Fx 

and 

h(x»y)  =  -(y+D  fxfxx  +  Fyy 

Then  one  has 

G(x,y)4»xx  +  Gx^x,y^x  +  ^yy  +  h(x>y)  =  0  (1) 

In  Reference  5  the  term  containinq  <fx  has  been  omitted. 

This  changes  the  details  of  the  argument,  but  not  the  overall 
conclusions.  Now  a  semidiscretization  is  carried  out  with 
respect  to  the  y  direction,  either  by  a  finite  element 
procedure  or  by  a  finite  difference  approximation  and  one 
arrives  at  a  system 


A(x)i|>"  +  A'fx)^'  +  +  f(x)  =  0  (2) 

Primes  denote  differentiation  with  respect  to  the  independent 
variable  x;  ip  and  f  are  x  dependent  vectors.  The  components 
of  are  the  parameters  which  serve  to  describe  the  potential 
along  lines  x  =  const.  (If  one  chooses,  for  instance, 
to  discretize  by  means  of  a  finite  difference  procedure. 


then  the  components  of  ^ (x)  are  given  by  the  values  of  |(x,y) 
for  the  pivotal  lines  y  =  const  chosen  for  the  discretization. 
In  this  case  the  components  of  f  are  the  values  of  h(x,y) 
along  those  lines.)  If  the  potential  along  a  line  x  =  const 
is  described  by  n  parameters,  then  xp  and  f  have  n  components. 

A  and  B  are  x  dependent  matrices.  (in  the  case  of  a  finite 
difference  approximation  A  is  a  diagonal  matrix  with  diagonal 
element  G(x,yn).  B  is  a  tri-diagonal  matrix  (with  elements 
l/h  ,  -  2/h  and  1/h  ) ,  it  arises  from  the  finite 
difference  form  of  <f> 

yy 

Incidentally,  Eq.  (1)  can  be  obtained  from  a  variational 
formulation,  namely 


$//( G(x,yH  2  +  *5  6  2  -  h(x,y)<(>)  dxdy  =  0 
x  y 


(with  suitable  boundary  conditions) .  The  differential  equation 
from  which  one  starts  in  Ref.  5  does  not  possess  this 
property.  The  system  Eq.  (2)  is  now  written  as  a  first  order 
system. 


Aip'  +  B\p  +  f  =  0 


where 


(3) 


(4) 


The  system  Eq.  (3)  can  obviously  be  solved  for  ip1  if  the 
determinant  | A |  =  | A |  is  different  from  zero.  Sinqular 
points  are  encountered  if  it  vanishes.  To  study  this  case 
consider  the  eigenvalue  problem 

A  u  -  Xu  =  0  (5) 


or  in  more  detail  with  u  = 
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U1  ~  ^U1  =  ® 


AU2  -  XU2  =  0 


It  has  n  trivial  eigensolutions 

X  =  +1,  U2  =  0,  u^  arbitrary 
Another  set  of  eigensolutions  has 


while  X  and  satisfy 

AU2  -  AU2  -  0 

The  further  development  is  best  carried  out  on  the  basis  of 
the  original  eigenvalue  problem,  Eq.  (5) .  We  characterize 
the  eigenvectors  and  eigenfunctions  by  subscripts  and  consider 
in  addition  the  adjoint  equation 


flun  -  Vr>  '  0 


*  -  Vi  ■ 0 

In  this  notation,  the  vectors  are  regarded  as  n  by  1  matrices, 
All  quantities  occurring  in  these  equations  depend  upon  the 
independent  variable  x.  The  determinant  of  the  matrix  A 
vanishes  at  those  values  of  x  where  one  of  the  eigenvalues 

vanishes.  Let  the  values  of  x  where  the  eigenvalue  A.  (x) 

k  K 

vanishes  be  denoted  by  x  .  One  has  (for  all  values  of  x) 

the  orthonormality  conditions 

’VI  “k  •  {»k  (8) 


Then  it  follows  that 


*1  A"k  =  Vilk 


We  represent  ip’  in  the  form 


ms 


(11) 


Then  one  finds,  using  Eqs.  (8)  and  (9) 

vj  (B$  +  f) 
k  Xk 

At  a  point  x  =  x1  ,  where  vanishes,  one  obtains  a  finite 

m 

value  of  1J7*  only  if 

*1  <»*)  =  0  (12) 

The  value  of  Ip  (xm)  is  determined  provided  that  Eq.  (12) 
which  restrict  the  choice  of  \p  is  satisfied,  although  Eq.  (11) 
gives  an  undetermined  expression  for  cm(xm).  Whether  or  not 
this  value  of  Ip  is  needed  depends  upon  the  method  used  for 
solving  the  systems  Eq.  (2)  or  Eq.  (3).  In  the  usual  predictor 
corrector  methods  one  evaluates  the  derivatives  at  the  chosen 
grid  points,  then  Ip  must  be  evaluated,  in  a  finite  element 
approach  the  derivatives  do  not  appear  explicitly,  and  one 
can  forgo  the  computation. 

The  derivative  I p  at  a  point  x  =  xm  can  be  computed  in 
the  following  manner.  We  define  the  auxiliary  quantity 

~i_  e  j? 

*  ■  xk 

Then  from  Eq.  (10) 


ip  =  ip  +  cu 
r  mm 


(14) 


One  has,  because  of  Eq.  (8) 


vT  tp’  =  0 
m  r 


For  x  =  xm  one  can  compute  ip  by  solving  the  equation 


(A+cuyJ)ip  =  -(Bipf) 

ip  =  -(A^u^)''  (B^+f) 


(15) 


F 


provided  that  Eq.  (12)  is  satisfied.  That  the  expression 

Eq.  (15)  is  correct  is  seen  because  the  matrix  (A  +  au  v  +) 

_  mm 

has  the  same  eigenvectors  as  A  and  also  the  same  eigenvalues 

except  for  k  =  m;  there  one  has  the  eigenvalue  X  +  a.  The 

m 

right  hand  side  of  Eq.  (15)  is  therefore  given  by 

r  vl  (fty+f)  J  (§*♦?) 

^  k^m  X^  uk  "  Xm+a  um 
The  second  term  vanishes  because  of  Eq.  (12) .  Now  we  must 
determine  cm  in  Eq.  (14).  One  obtains,  by  applying  L'Hospital's 
rule  to  Eq.  (11) 


c 


m 


(dV|J  /dxHify+f)  +  vj  ( (dB/dx)^  +  (df/dx))  +  B(diJ)/dx) 

dxydx 


(16) 


To  compute  the  quantities  occurring  in  this  expres 
differentiate  Eq.  (7)  : 


dvT  _ 
m  A 


dvT 

m 


X  +  vT 
m  m 


^  -  p"  vT  =0 

dx  dx  m 


dx  n  dx  . 

The  solvability  condition  for  this  inhomogeneous  system  is 


Hence 


(dA  .  dXm,  - 
m  vdx  dx  '  m 


dX  -  aI 

_ m  _  T  dA  - 

dx  vm  dx  um 


(17) 


Furthermore , 


since 


■v  I  m\ 

X  (x  ) 
m 

dvj  - 
m 


dx 


A  = 


0 


vT 

m 


(M  d\ru 
ldx  ‘  dx  ' 


Hence,  employing  an  argument  similar  to  the  one  used  above 

dX 


dvT 

m 


dx 


=  .  ;t  (dA 
vm  vdx 


dx 


— ) 


(A+au  vj)‘ 
m  m 


+  cv. 


(18) 


Here  c  is  arbitrary  (ultimately  determined  by  a  normalization 
condition) .  The  constant  c  vanishes  upon  substitution  of 
Eq.  (18)  into  Eq.  (16),  because  of  Eq.  (12).  Multiplying 
Eq.  (16)  by  dXm/dx  and  substituting  Eq.  (14)  one  finds 
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c™  C5m  3?  “ml  pm  (S7  ■  (5m  f>.  °m,)(A  *  “Vi  >“<5*fh 

t5J  (f*  *  -  5J  «***%,»;  *  Vl  BSm  • 0 


Hence 


cm  •  C5T  c  *  B>5J"'  <*I  <£  -  <vT  £  =m> *  B)«^v.  )"(5*t'f) 


-vT  (^  +#)} 
m  vdxv  dx' 

This  formula  expresses  c  in  terms  of  and  the  known 

m 

quantities  A,  B,  dA/dx,  dB/dx,  f  and  df/dx.  As  was  mentioned 
above,  ip  is  not  always  needed. 

To  impose  the  condition  that  the  solution  of  the  system 
is  free  of  singular  points,  one  must  first  determine  the 
values  of  x  for  which  =  0.  For  these  points  one  must 
evaluate  u^  and  v^  and,  for  methods  in  which  ^ '  is  needed,  a 
number  of  derivatives  of  known  quantities. 

For  an  illustration,  consider  the  second  order  equation 


With 


-x<{>"  -<J>'  -<}>  +  f  =  0 


i^l  =  <t>.  4>2  =  4* 


one  obtains  the  equivalent  first  order  system 


i  o  rv 


0  -xi  uy 


o  -l  U 


i  * 


The  associated  eigenvalue  problem,  Eq.  (5) ,  is  given  by 


The  eigenvalues  and  eigenfunctions  are  found  by  inspection 


x1  *  1; 

S' 

V 

» 

p  i 

vi 

* 

« 

i 

•u2* 

1 

.0. 

-V2- 

i 

:o. 

^2  *  “ x  * 

S1 

- 

V 

* 

w 

i 

o" 

•u2- 

O 

J  - 

U2. 

i 

ul- 

A  singular  point  occurs  at  x  =  0,  (A2(x)  is  zero  at  this 

point).  The  regularity  condition  Eq.  (12)  gives 


or 


=  0 


[0  :  l]  f-*2 

-<J>2  +  f  J 

-1^(0)  -  4>2( 0)  +  1(0)  =  0 

-<j>( 0)  -  4>‘(0)  +  f(0)  =  0 


(21) 


This  result  can  be  found  directly  from  Eq.  (19)  .  The 
derivative  ip^'  (0)  computed  from  Eq.  (20), 

-x  ^2’  -  ^2  "  'h  +  f  =  ® 

is  undetermined.  The  desired  result  is  obtained  by  a 
differentiation  with  respect  to  x 

i\>2  (0)  =  (l/2)(-*1'(0)  +  f'  (0)) 
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moreover 


V(0)  *  ♦2(0) 

The  essence  of  a  numerical  procedure  which  computes  the 
solution  in  an  interval  -a  <  x  <  0,  with  4>  ( —a )  prescribed 
is  already  seen  if  one  traverses  this  interval  in  one  step. 

One  has  in  the  simplest  difference  formulation  of  Eq.  (19) 

a  j>(-a)-2(j>(-a/2)+<|>(P)_  _  <fr(0H(-al  .  ^(-a/2)  +  f (-a/2)  =  0 
2  (a/2)2  a  (23) 

The  value  of  <M-a)  is  given  by  the  boundary  conditions.  In 
addition  one  has  the  regularity  condition  for  x  =  0. 

It  is  consistent,  if  one  expresses  the  value  of  <t'(0)  by  a 
finite  difference  approximation. 

4.(0)  =  34>(0)-4<t>(-a/2)-Hl>(-al 

9  1 '  a  (24) 

The  values  of  <H-a/2)  ,  0(0),  and  <p’(0)  can  then  be  computed 
from  Eqs.  (23),  (24)  and  (21). 

In  the  first  step  of  the  continuation  to  positive  values 
of  x  a  formula  corresponding  to  Eq.  (24)  is  used 

,(0)  =  -34>(0)+4<}>(a/2)-4>(a) 

d 

The  left  side  is  given  by  Eq.  (24).  In  addition,  one  has 
from  the  difference  approximation  of  the  differential  equation 
for  the  point  x  =  a/2 

_a  j>(0)-2<l>(a/2)+(l>(d)  _  4>(a)-<t>(0l  _  ^(a/2)+f (a/2)  =  0 
2  (a/2)2  a 

The  last  two  equations  are  used  to  compute  <Ma/2)  and  0(a). 

From  there  on  the  values  of  0  can  be  computed  in  sequence. 

The  procedure  is  even  simpler  if  the  grid  points  straddle 
the  point  x  =  0.  Assume  that  0 ( —  3/4a)  is  prescribed.  Then 
one  has,  for  the  first  inner  point 

12 


**Mr  ^  H3PW 


^  - 


-  4>( -a/4)+f ( -a/4)  *  0 


a  <t(3a/4)-2<t>(-a/4)+<t>(a/4)  _  6(a/4)-<fr(-3a/4) 

1  un?  * 

In  the  regularity  condition,  derived  from  Eq.  (19)  ,  <J)q  is 
expressed  either  by  linear  interpolation 

<J>(0)  =  (l/2)[<J>(-a/4)+<fr(+a/4)] 

or  by  quadratic  interpolation 

4>(0)  =  -( l/8)<J>(-3a/4)+(3/8)^(a/4)+(3/4)<|>(-a/4) 

Moreover 

4>  *  (0 )  =  ^a/-4^ra/'-^ 

One  thus  has  two  equation  which  allow  one  to  compute 
4>(-a/4)  and  <M+a/4)  .  For  a  continuation  to  further  positive 
values  of  x,  one  uses  a  marching  procedure.  The  value  of 
( 3a/4 )  is,  for  instance,  obtained  in  terms  of  <j>(-a/4)  and 
<M+a/4)  ,  from  the  difference  equation  for  the  point  x  =  a/4. 

The  situation  is  similar  if  one  treats  Eq.  (20)  instead 
of  Eq.  (19)  by  a  difference  formulation.  Taking  as  grid  points 
x  =  -a,  x  =  -a/2  and  x  =  0,  one  obtains  as  equations  for  the 
middle  of  the  intervals  -a  <  x  <  -a/2  and  -a/2  <  x  <  0 


^1  (-a) 

1/2 


t^2(-a/2)+^2(-a) 

2 


=  0 


1^(0)-^  (-a/2)  <|>2(0)+*2(-a/2) 

- a'/2 - 2  =  0 

3a  (-a/2)-H/>1  (-a) 

4  a72  2 


ii»2(-a/2)+i(i2(-a) 

2 


=  0 


a  4>2(0)-g;2(-a/2)  ^1(0)+^1(a/2)  ^2(0)+^2(-a/2)  / 

? - i72 - 2 - 2  +  f(‘  V  *  0 


The  value  of  t^^(-a)  is  given  as  boundary  condition.  In 
addition,  one  has  the  regularity  condition  (first  of  Eqs.  (21)). 
One  thus  has  6  equations  for  6  unknowns  (ij>.  and  eac^ 
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of  the  points  x  =  -a,  x  =  -a/2  and  x  =  0) .  For  the  continuation 
to  the  right  one  uses  a  marching  procedure  which  considers 
one  interval  at  a  time. 

No  essential  differences  are  encountered,  if  one  obtains 
discretized  equations  by  a  finite  element  procedure.  If 
one  takes,  for  instance,  third  degree  shape  functions,  then 
one  characterizes  the  solution  by  the  parameters  <p(-a), 

(-a),  <J>(0)  ,  <t>'(0),  etc.  For  each  interval  one  has  two 
weight  functions  (because  one  has  tv/o  parameters  for  each 
grid  point).  In  addition,  <p(-a)  is  assigned  and  one  has 
the  regularity  condition  (second  of  Eqs.  (21)).  This  suffices 
to  compute  the  solution  for  x  <  0. 

Eq.  (22)  for  the  derivative  at  x  =  0  would  be  needed 
if  one  uses  formulae  analogous  to  those  occurring  in  the  usual 
integration  techniques  for  ordinary  differential  equations. 

Using  a  very  simply  integration  procedure  and  traversing  the 
interval  -a  <  x  <  0  in  two  steps  one  has 

*2(-a/2)-iP2(-a)  -  (a/2)(l/2)(ijy(-a/2)+V(-a)) 

4»1(-a/2M1(-a)  *  (a/2)(l/2)(^  ’  (-a/2)+^ '  (-a)) 

U>2(0)  -  u»2(-a/2)  -  (a/2)(l/2)(^2'  (0)+^2*  (-a/2)) 

^(0)  -^(-a/2)  =  (a/2) (1/2) (^-i' (0)+^'  (-a/2)) 

At  the  points  -a  and  -a/2,  ip^'  and  ipj*  are  expressed  in  terms 
of  ip^  and  4>2  by  means  of  the  differential  equations  (20).  At 
x  =  0,  these  quantities  are  expressed  by  Eqs.  (22).  They  are 
derived  from  Eq.  (20).  Thus,  one  has  four  relations  for  six 
quantities  (ip^  and  ip2  at  three  points),  tp^(-a)  is  given  by  the 
boundary  condition.  In  addition,  one  has  the  regularity 
condition  (21). 

Next,  we  consider  an  example  which  arises  from  a  partial 
differential  equation.  Consider 
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(25) 


(y-a(x))<j>xx-a'(x)<j>x+4>yy  =  0 


For  y  >  a(x)  the  problem  is  elliptic,  for  y  <  a(x)  it  is 
hyperbolic.  The  line  y  =  a(x)  will  be  called  the  parabolic  line. 
We  choose  as  boundary  condition,  4>  =  0  for  y  =  0,  and  y  =  1. 

Let  a(x)  be  a  function  which  changes  in  the  region  under 
consideration  monotonically  from  some  value  less  than  zero 
to  a  value  exceeding  1.  Then  one  has  values  of  x  at  the  left 
where  the  problem  is  elliptic  and  at  the  right  where  it  is 
hyperbolic.  In  the  intermediate  region  one  finds  a  parabolic 
line  which  goes  from  the  lower  left  to  the  upper  right.  Let 


y  =  - 

N 


(n  <  N)  .  We  approximate  <p(x,y)  by  means  of  the  values  which 


it  assumes  along  the  "pivotal"  lines,  y  =  y  . 


4>n(x)  =  4>(x,yn) 


Taking  a  finite  difference  approximation  for  <|>  ,  one  then 

obtains  the  system  of  ordinary  differential  equations 


(yn-a(x))*;  -  a'(x)^  +  N~2(<t>n+i (x)  -  2<t>n(x)  +  4>n_-| (x) )  =  0 


n  =  1,  2  ...  N  -  1 


The  matrix  A  introduced  in  Eq.  (1)  is  a  diagonal  matrix  with 


elements  y^  -  a(x),  and  the  augmented  matrix  A  is  also  a 


diagonal  matrix  with  additional  elements  1.  Eigenvalues  0 


occur,  obviously,  at  values  of  x,  to  be  denoted  by  x,  where 


a(xn)  s  yr 


These  are  the  intersection  of  the  pivotal  lines  with  the 
parabolic  line.  It  is  obvious  from  the  form  of  the  differential 
equations,  that  these  are  the  singular  points.  Let 


r'*t"  -pm  ^ 


i 


I; 


l 


nth  component 


be  the  unit  vector  in  the  direction  of  the  n^*1  component  of 
the  space  of  the  vectors  <J>.  One  has  as  eigenvectors  of  A»  aside 
from  those  with  eioenvalue  1  which  are  not  of  interest 


0 

0 

u  = 

e 

m 

»  vm  = 

*  m 

*3* 

e 

m 

L  J 

with  eigenvalues  X  (x)  =  y  -  a(x) . 

m  m 

the  matrix  B 


In  this  case  one  has  for 


B  =  N 


-2 


■2  1 
1  -2  1 


2  1 
1  -2 


The  regularity  condition  (12)  at  a  singular  point  x 
gives 


=  -a' 


m 


=  x 


1  CO 

i  a? 

o 

*1 

o 

ii 

-♦2 

^2 

L  J 

~a'  (xm)\p2 '  +  B^i 

Hence 


-a'(xmH2,m  +  rtW**)  "  2Wxm>  +  Vm+l(xm)]  =  0  (27) 


This  condition  can  also  be  found  by  inspection  of  Eq.  (26) . 

One  needs  2n  boundary  conditions  for  the  system  of  ordinary 
differential  equations  Eq.  (26)  .  n  conditions  are  given  by 
the  values  of  <J>(x£,y  ),  where  x^  is  the  value  of  x  at  the 
left  end  of  the  region  under  consideration,  n  further 
conditions  are  obtained  from  the  regularity  conditions  Eq.  < 27). 
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The  form  of  the  regularity  conditions  is  more  complicated 
if  one  applies  a  different  form  of  discretization,  although 
it  is  based  on  the  same  idea.  One  might  choose  a  finite 
element  approximation  for  the  y  direction,  or  (even  though 
this  may  not  be  practical)  a  representation  in  terms  of  a 
finite  Fourier  analysis  based  on  the  values  of  <J>  (see  Ref.  5) 
In  those  cases,  the  matrix  A  does  not  have  diagonal  form,  and 
the  determination  of  the  eigenvalues  and  eigenvectors  vn  is 
more  difficult.  In  Ref.  5  it  is  shown  by  an  example  that  for 
large  values  of  N  (that  is  if  the  pivotal  lines  y  =  const 
lie  close  together)  the  functions  of  y  represented  by  the 
eigenvectors  of  A  approach  6  functions.  Thus,  one  obtains  in 
essence  the  same  eigenfunctions  as  in  a  finite  difference 
discretization.  They  are  used  again  to  derive  regularity 
conditions  for  the  points  x  =  xn. 

The  system  of  ordinary  differential  equations  so 
obtained  is  now  treated  by  a  discretization  with  respect 
to  x  by  the  same  techniques  described  above  for  a  simpler 
case. 

The  weight  functions  applied  in  the  x  discretization 
combined  with  the  weight  functions  for  the  y  discretization 
give  two  dimensional  weight  functions.  Figures  1  and  2 
show  the  boundaries  of  the  regions  to  which  the  weights  are 
applied.  For  constant  weights  the  boundaries  of  adjacent 
areas  coincide.  If  one  works  with  nonconstant  weights  then 
there  is  in  general  an  overlap  of  the  weight  areas,  but  the 
general  arrangement  is  the  same. 

Figure  1  refers  to  bilinear  shape  functions  which  are 
characterized  by  the  values  of  <J>  at  the  corners  of  the 
quadrangular  elements.  To  each  of  the  inner  lines  y  =  const 
which  form  an  element  boundary  there  belongs  one  differential 
equation  and  for  each  of  these  lines  one  needs  one  regularity 
condition.  In  the  supersonic  region  it  is  probably  necessary 
for  reasons  of  stability  to  shift  the  weight  areas  downstream. 
This  has  not  been  indicated  in  the  figure.  Each  row  of 
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Figure  1.  Weight  Areas  for  Bi-Linear 
Shape  Functions,  Element 
Boundaries  and  Sonic  Line. 
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Figure  2.  Weight  Areas  for  Bi-Cubic  Shape 
Functions,  Element  Boundaries 
and  Sonic  Line. 
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weight  areas  correspond  to  one  differential  equation.  (The 
differential  equations  will,  of  course,  contain  terms  which 
originate  from  adjacent  lines  y  =  const.)  For  each  of  the 
unknowns,  that  is  for  each  row  of  weight  areas  one  has  one 
regularity  condition.  For  Eq.  (1)  it  is  the  requirement 
that  Gx<J>x  +  4>yy+  h(x,y)=  0  at  the  parabolic  line. 

For  bicubic  elements  the  same  situation  is  shown  in 

Fig.  2.  The  lengths  of  the  sides  of  the  quadrangles  is 

taken  twice  those  of  the  bilinear  elements.  Bicubic  elements 

can  be  characterized  by  the  values  of  <J>,  $  ,  6,  and  $  at 

x  y  xy 

each  element  corner.  For  each  inner  corner  one  therefore 

has  four  weight  areas.  At  the  upper  and  lower  boundaries 

the  values  of  4>  are  known,  in  addition  one  knows  at  the  left 

boundary  the  values  of  and  at  the  upper  and  lower  boundaries 

the  values  of  <{>  ,  One  has  therefore  two  weiqht  areas  for  each 

element  corner  at  the  boundary  of  the  reqion,  Aqain,  each 

row  of  weight  areas  corresponds  to  one  differential  equation. 

And  for  each  such  row  one  obtains  from  the  system  of 

ordinary  differential  equation  one  regularity  condition. 

In  the  fotm  derived  above  these  conditions  are  inconvenient. 

They  require  that  one  determine  the  values  of  x11  (those  values 

of  x  where  the  eigenvalues  A (x)  vanish) ,  and  the  eigenfunction 

vn  of  the  adjoint  operator  A.  The  exact  form  of  these  conditions 

depends  upon  the  chosen  discretization  procedure.  These 

conditions  express  in  discretized  form  the  fact  that  (for 

the  present  example)  — G  <J>  +  <J>  +  h(x,y)  =  0  at  the  parabolic 

x  x  yy 

line.  The  procedure  derived  above  is  one  possible  form  which  it 
can  be  expressed.  It  is  simpler  to  introduce  this  requirement 
directly.  Details  of  such  a  procedure  are  shown  in  the  next 
section. 

If  one  wants  to  insure  overall  conservation  of  mass, 
then  one  must  use  constants  as  weight  functions  and  cover  the 
entire  region  of  integration  with  weight  areas,  neither  gaps 
nor  overlap  are  permitted.  Usually,  conservation  of  mass  is 
needed  only  for  shocks. 
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In  Ref.  5  the  question  has  been  raised  of  what  happens 
if  lines  x  =  const  intersect  certain  characteristics  more 
than  once.  For  a  more  complicated  partial  differential 
equation  this  can  easily  happen  in  the  vicinity  of  the  parabolic 
line,  where  the  two  families  of  characteristics  have  a 
common  tangent.  If  one  uses  in  such  a  case  the  semidiscretiza¬ 
tion  technique  described  above,  then  one  obtains,  again,  a 
system  of  differential  equations  with  singular  points,  but 
these  singular  points  lie  in  the  supersonic  region.  The 
typical  phenomena  associated  with  the  parabolic  line  are 
postponed  by  some  kind  of  a  sweep  effect  introduced  by  the 
orientation  of  the  coordinate  system.  This  approach  is 
questionable  because  one  can  find  examples  in  which  it  seems 
to  give  solutions  to  boundary  value  problems  which  are  not 
well  posed. 

The  question  whether  it  is  necessary  to  choose  the 
element  boundaries  so  that  this  phenomenon  cannot  occur,  or 
whether  it  suffices  if  one  chooses  the  boundaries  of  the 
weight  areas  so  that  the  regions  of  dependence  are  properly 
taken  into  account  in  the  weighting  of  the  residuals  is  left 
open. 
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SECTION  III 


DIRECT  TREATMENT  OF  THE  CONDITIONS 
FOR  THE  SONIC  LINE 

The  discussion  in  Section  II  started  from  the  differential 
equation  (1)  obtained  by  linearization  of  the  transonic  version 
of  the  potential  equation 

G(x,y)<J>xx  +  Gx(x,y)<t>x  +  <t>yy  +  h(h,y)  =  0 

A  discretization  is  carried  out  by  introducing  a  rectangular 
mesh,  expressing  <p  by  either  bilinear  or  bicubic  shape 
functions  and  by  applying  weights  (say  constant  weights) 
according  to  the  patterns  shown  in  Figures  1  and  2,  respectively 
for  bilinear  and  bicubic  shape  functions.  Along  the  parabolic 
line,  that  is  along  the  line  where  G(x,y)  =  0,  the  differential 
equation  yields  the  relation 

G  <j)  +  <{>  +  h(x,y)  =  0  (28) 

x^x  yy 

Eq.  (28)  will  be  used  to  derive  further  equations  in  addition 
to  those  obtained  by  weighting  of  the  residual.  One  can  assume 
that  the  weight  areas  shown  in  Figures  1  or  2  arise  in  two 
steps  by  carrying  out  a  discretization  in  the  y-direction  first 
and  in  the  x-direction  afterwards.  In  the  first  step  one 
obtains  a  system  of  ordinary  differential  equations,  one  for 
each  row  of  weight  areas.  For  each  differential  equation  and 
therefore  also  for  each  row  of  weight  areas  one  needs  one 
additional  condition.  These  are  obtained  from  Eq.  (28) .  These 
conditions  are  written  in  the  form 

(29) 

/(Vx  +  ‘t’yy  +  h(x.y))dS  =  0 

where  the  integrals  are  formed  over  the  segments  of  the 
parabolic  line  which  extend  over  the  width  of  the  rows  of 
weight  areas.  This  provides  as  many  conditions  as  there  are 
differential  equations.  The  discussions  of  Section  II  have 
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the  purpose  of  showing  that  these  conditions  arise  naturally 

from  the  system  of  differential  equations  obtained  by  the 

y-discretization.  Proceeding  in  this  manner,  one  obtains 

the  system  symbolized  for  bilinear  shape  functions  by 

Fig.  3.  The  shape  functions  are  determined  by  the  values  of 

$  at  the  corners  of  the  quadrangles.  These  values  have  been 

denoted  by  <J> .  ,  where  j  is  the  subscript  for  the  element 
3  t K- 

boundaries  x  =  const,  and  k  for  the  element  boundaries  y  = 

const.  The  values  of  ,  .  for  constant  j  (that  is  the  values 
th  J ,K 

of  <t>  along  the  j  line  x  =  const)  are  combined  into  a  vector 
ipj .  Part  of  the  matrix  shown  in  Fig.  3  is  block  tridiagonal, 
the  blocks  are  square  with  dimensions  determined  by  the 


dimension  of  the  vector  \p ^ .  The  first  row  contains  only  two 
blocks  because  at  the  entrance  cross  section  the  vector  ip  is 


given  by  the  boundary  conditions,  its  contribution  to  the 
system  of  equations  appears  in  the  inhomogeneous  part.  The 
last  row  of  the  square  blocks  contains  three  because  the 
vector  ip  in  the  exit  cross  section  is  unknown.  The  last  row 
consists  of  single  equations.  They  arise  from  the  regularity 
conditions  Eq.  (28).  The  segment  of  the  sonic  line  which  lies 
within  one  row  of  weight  areas  may  intersect  more  than  one 
weight  area.  Accordingly,  the  conditions  expressed  by  Eq. 

(28)  will  give  relations  between  a  number  of  vectors  t|/j  ,  but 
usually  only  a  very  small  number  of  components  of  these 


vectors  will  appear. 


The  sparsity  of  this  matrix  can  be  taken  into  account 
fairly  easily.  In  principle,  one  eliminates  the  vectors  ^ 
in  the  sequence  given  by  their  subscripts.  The  elimination 
includes,  of  course,  the  equations  obtained  from  the 
regularity  conditions  Eq.  (28) .  In  this  process  there  comes  a 
points  where  some  equations  obtained  from  the  regularity 
conditions  contain  after  the  eliminations  only  components  of 
the  first  of  the  remaining  vectors  iji ^ .  It  is  then  possible 
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to  express  some  of  the  components  of  this  vector  in  terms  of 
the  other  components.  While  in  the  initial  elimination  process 
one  has  to  solve,  in  each  step,  systems  of  equations  whose  size 
is  determined  by  the  dimension  of  ,  one  performs  from 
this  point  on  this  task  in  two  steps;  in  the  first  one,  some 
of  the  components  of  the  vector  iJk  are  expressed  in  terms 
of  the  others,  in  the  second  one,  one  eliminates  the  remaining 
components  from  the  system.  After  all  regularity  conditions 
have  been  taken  into  account  in  this  manner,  one  arrives  at 
a  system  in  which  one  can  compute  the  remaining  vectors 
directly  in  sequence.  This  means  that  one  carries  out  a 
marching  procedure. 
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SECTION  IV 


TRANSITION  FROM  A  SUPERSONIC  TO  A 
SUBSONIC  FLOW 


For  a  first  orientation  we  consider  ,  again,  the  ordinary 
differential  equation  (19) 


-x4>"  -<J>'  -<J>  =  0 


or  rather 


x<j>"  +  4>  ’  -<)>  =  0 


and  with  a  slight  generalization 

a(x)4>"  +  a'(x)<j>'  -<p  =  0,  a(0)  =  0,  a'(x)>0  for  x>0 


(30a) 


because  it  is  natural  to  make  the  transition  from  the 
oscillatory  type  of  the  particular  solutions  (which  is 
analogous  to  a  supersonic  flow)  to  the  nearly  monotonic  type 
(analogous  to  subsonic  flows)  as  one  proceeds  in  the  positive 
x  direction  from  a  negative  to  a  positive  value.  The  left 
end  of  the  region  in  which  the  solution  is  to  be  found  is 
taken  at  a  negative  value  of  x.  There  one  prescribes  the 
values  of  <j>  and  <}> '  ,  in  analogy  to  boundary  conditions  for 
supersonic  problems.  The  right  end  of  the  region  is  taken 
at  some  positive  value  of  x,  and  in  accordance  with  boundary 
conditions  for  elliptic  problems  one  prescribes  only  the  value 
of  <f>  (or  * )  .  While  we  had  too  few  boundary  conditions  in  the 
discussions  of  Section  II  (1  instead  of  2)  ,  we  have  now  too 
many  (3  instead  of  2).  The  differential  equation  has  a  singular 
point  at  x  =  0.  We  found  in  Section  II,  that  for  the  transition 
from  the  "elliptic"  to  the  "hyperbolic"  type  the  missing 
condition  is  provided  by  postulating  regularity  of  the  solution 
at  x  =  0.  Under  the  present  condition  regularity  cannot  be 
achieved  because  for  negative  values  of  x  the  solution  is 


completely  determined  by  the  initial  conditions  at  the  left 
end  of  the  interval.  We  ask  whether  one  can  satisfy  the 
boundary  conditions  at  the  right  end,  if  at  x  =  0  one  satisfies 
the  differential  equations  only  in  the  weak  sense. 

In  Eq.  (30a)  the  coefficient  a(x)  is  regular  and  has 
a  zero  of  the  first  order  at  x  =  0.  Particular  solutions  of 
Eq.  (30a)  can  be  written  in  the  following  manner 

=  P^x)  (31) 

and 

<P2  =  log]  x|  Pj(x)  +  cxP2(x)  (32) 

where  P^(x)  and  P2(x)  are  Power  series  in  x  which  have  1  as 
constant  term,  and  c  is  a  suitable  constant.  Incidentally, 

Eq.  (30)  is  solved  by 

<t>  =  x  >  0 

<f>  =  2Q( | xf*5)  x  <  0 

where  Zq  is  a  linear  combination  of  the  Bessel  functions 
JQ  and  Nq.  The  general  solution  is 

4>(x)  =  c-,+4>1(x)  +  c2+<t>2(x)  ,  x  >  0  (33) 

4>(x)  =  c-j tJ>-j ( x )  +  c2^2^ x ^  »  x  <  0 

cL+,  c2+,  c"  ,  c2'  are  constants.  In  the  region  x  <  0  the 
solution  and  therefore  also  c^-  and  c2~  are  determined  by 
the  conditions  prescribed  at  the  left  end  of  the  interval. 

Eq.  (30a)  can  be  written  as 

^  (a(x)<t>; (x))  -  4>  =  0 

Then  by  integrating  from  x  =  -e.  to  x  =  +e2  (e.  >  0,  e2  >  0) 


(34) 


a(e2)<f>' (e2)  -  a(-e1  )<J>' 


*  0 


The  integral  vanishes  in  the  limit  0,  e2  0,  if  one 

expresses  <p  by  Eqs.  (33)  ,  (31)  and  (32)  . 

To  evaluate  the  first  two  terms  we  form 

<J>'(x)  =  C-| ~P -| ' (x)  +  c2+  [1  P1  +  log  x  P1<  +  cP2  +  cxP2‘] 

Then,  since  a(0)  =  0,  and  since  P^(0)  =  1 

11m  (a(e)<f>( (e) )  =  cj1  a '  (0) 

e-*0 

Therefore,  from  Eq.  (34) 
c2  a  c2 

while  c^  remains  undetermined. 

The  requirement  that  at  x  =  0  the  differential  equation 
be  satisfied  only  in  the  weak  sense  introduces  an  indeterminacy 
which  makes  it  possible  to  satisfy  the  boundary  condition 
prescribed  at  the  right  end  of  the  interval. 

The  numerical  procedure  is  quite  clear  if  one  treats 
the  differential  equation  by  some  standard  numerical  approach 
for  the  solution  of  ordinary  differential  equations.  For 
negative  values  of  x,  one  has  an  initial  value  problem  which 
is  solved  up  to  some  negative  value  of  x  close  to  zero.  In 
the  vicinity  of  x  =  0  one  computes  two  linearly  independent 
particular  solutions  for  <j>  by  means  of  their  development 
and  determines  the  values  of  c^  and  c2  by  matching  at  some 
negative  value  of  x  with  the  solution  obtained  by  numerical 
integration.  For  positive  values  of  x,  one  starts  from  x  -  0, 
setting  c2+  =  c2~  and  leaving  c^+  arbitrary;  c^+  is  determined 
by  the  boundary  conditions  at  the  right  end. 

In  studying  which  form  the  procedure  assumes  if  the 
differential  equations  are  solved  by  a  finite  element  approach. 
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let  us  assume  that  the  elements  straddle  the  point  x  =  0. 
Choosing  the  weight  regions  for  x  <  0  in  a  manner  suitable 
for  hyperbolic  problems  and  for  x  >  0  in  a  manner  suitable 
for  elliptic  problems,  the  number  of  equations  exceeds  the 
number  of  unknowns  by  one  (Fig.  4a) .  The  weighting  of 
residuals  is  a  finite  dimensional  approximation  to  the  concept 
of  weak  solutions,  the  present  procedure  is  basically  in 
conformance  with  the  concept  developed  above,  but  it  gives 
too  many  conditions.  This  is  remedied  by  omitting  the  weight 
region  which  covers  the  first  finite  element  boundary  at  a 
positive  value  of  x  and  by  extending  the  preceding  weight 
area  so  that  it  covers  this  boundary  (Fig.  4b) .  This  treat¬ 
ment  of  the  vicinity  of  x  =  0  is  the  analogue  to  Murman's 
shock  operator. 

The  linearized  form  of  the  transonic  equation  Eq.  (1) 
which  served  as  basis  for  the  discussions  carried  out  so  far, 
is  complicated  because  the  function  G  in  Eq.  (1)  is  discon¬ 
tinuous  at  the  shock  location.  For  a  first  orientation  we 
therefore  consider  the  nonlinearized  transonic  equation 

-(y+D^xX  +  V  =  0 

Fig.  5  show  the  choice  the  weight  areas  for  such  a  case. 

For  bilinear  elements  one  has  delta  function  residuals  at 
the  element  boundaries.  The  weight  areas  must  overlap  the 
lines  at  which  delta  functions  occur.  In  the  supersonic  region 
one  best  uses  an  infinitesimal  overlap  at  the  upstream  end 
of  the  elements,  but  not  at  the  downstream  end.  In  the 
subsonic  region,  the  weight  areas  are  chosen  symmetric  (in 
the  x  and  in  the  y  directions)  with  respect  to  the  element 
boundaries.  In  elements  which  contain  the  shock,  the  weight 
areas  overlap  the  downstream  as  well  as  the  upstream 
boundary.  The  first  weight  areas  of  subsonic  type  start  at 
the  downstream  boundary  of  the  weight  area  pertaining  to  the 
shock.  In  the  shock  weight  areas  the  weight  is  taken  constant 
to  guarantee  conservation  of  mass. 
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Figure  4a.  Weight  Areas  for  an  Ordinary  Differential 
Chosen  According  to  the  Type  (Oscillatory 
or  Nonoscillatory)  of  Particular  Solutions 
(More  Equations  than  Unknowns) . 
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Figure  4b.  Modification  of  Figure  4a  in  which  one 
Weight  Area  is  Eliminated. 


Figure  5.  Weight  Areas  for  the  Transition  From 
Supersonic  to  Subsonic  Velocities  by 
Means  of  a  Shock. 


In  examining  whether  mass  is  conserved  through  the 
shock,  one  is  hampered  by  the  fact  that  shape  functions 
which  are  cutomarily  applied  are  not  suited  to  express  the 
sudden  change  of  velocity  which  occurs  in  a  shock.  Conservation 
of  mass  can  be  shown  only  in  the  average  over  several  rows 
of  weight  areas,  and  this  requires  that  the  states  in  front 
and  behind  the  shock  change  only  little  along  the  shock  and 
that  the  location  of  the  shock  relative  to  the  weight  areas 
in  the  upper  and  in  the  lower  rows  is  the  same.  Then  the 
inadequacies  in  the  representation  of  the  flow  field,  in 
particular  of  <J>  ,  cancel  each  other  and  one  obtains  that, 
in  the  average,  mass  is  conserved.  This  is  exactly  Murman's 
argument . 

The  shock  operator  is,  of  course,  compatible  with  a 
smooth  transition  from  a  supersonic  to  a  subsonic  flow.  It 
is  therefore  applicable  to  shockfree  flow  fielus. 

The  fact  that  the  shock  conditions  are  satisfied  only 
in  the  average  over  several  elements  makes  the  choice  of 
rather  small  elements  mandatory.  At  least  in  the  staue  of 
the  computations  where  a  relatively  good  approximation  has 
already  been  found  and  one  wants  to  refine  this 
by  using  elements  in  which  the  shock  position  is  identified. 

A  shock  is  without  influence  on  the  supersonic  region  upstream. 
The  supersonic  region  can  actually  be  continued  downstream 
beyond  the  shock  location.  (An  obvious  example  is  the  bow 
shock  in  front  of  a  blunt  body  in  a  supersonic  parallel 
flow.)  In  the  supersonic  flow,  the  state  of  points  which 
lie  in  the  area  of  influence  of  the  data  at  a  preceding  station 
x  =  const  is  independent  of  the  boundary  conditions  which  one 
may  prescribe  at  the  upper  or  lower  boundary  of  the  region. 

This  is  approximately  correct,  even  if  one  applies  it  in  an 
implicit  method,  as  is  frequently  done  in  this  context. 

One  therefore  can  find  the  supersonic  reqion  ahead  of  the 
shock  and  to  some  extent,  its  continuation  downstream  with 
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the  usual  finite  element  formulation  for  supersonic  flows 
using  fictitious  boundary  conditions  at  the  upper  end  (this 
description  refers  to  a  shock  on  the  upper  side  of  a  profile) . 
Once  the  flow  in  the  supersonic  region  ahead  of  the  shock  is 
known,  the  complete  state  of  the  flow  immediately  behind  the 
shock  depends  only  upon  the  position  of  the  shock.  For  the 
present  discussion,  we  are  using  rectangular  elements  and 
consider  in  lieu  of  ordinary  differential  equations  the 
difference  equations  corresponding  to  rows  of  weight  areas 
bounded  by  lines  y  =  const.  The  shock  can  be  described  by 
means  of  the  x  coordinates  of  its  intersection  with  these 
row  boundaries.  These  coordinates  are  now  considered  as 
unknown.  In  the  simplest  case  one  can  approximate  the  shock 
between  the  intersection  points  by  straight  lines.  In  this 
manner  one  obtains  one  additional  free  parameter  for  each  row 
of  weight  areas.  This  gives  exactly  the  number  of  additional 
unknowns  to  satisfy  the  boundary  conditions  at  the  downstream 
end  of  the  subsonic  region. 

An  insight  into  the  character  of  the  problem  can  be 
obtained  in  the  following  manner:  Assume  that  the  shock 
position  is  known.  Then  one  has  in  the  region  downstream  of 
the  shock  a  purely  subsonic  problem.  A  finite  element 
formulation  will  lead  to  a  system  of  equations  which  contains 
among  the  unknowns  the  potential  at  the  shock,  while  the  mass 
flux  through  the  shock  (which  is  temporarily  assumed  to  be 
known)  appears  as  a  contribution  to  the  inhomogeneous  terms. 
The  solution  to  this  problem  would  be  completely  determined. 
Actually,  the  normal  component  to  the  mass  flow  is  unknown 
because  the  shock  location  is  not  known,  but  it  and  also  the 
potential  can  be  expressed  in  terms  of  the  shock  abscissas. 

The  potential  is  continuous  and  normal  component  of  the  mass 
flux  are  continuous  through  the  shock.  This  qives  two  sets 
of  additional  conditions  in  terms  of  the  shock  abscissas 
which  are  considered  as  additional  parameters  in  the 


description  of  the  flow  field.  The  procedure  is  best  used 
in  the  form  of  a  Newton  Raphson  iteration  for  the  whole  flow 
field  in  which  one  computes  simultaneously  the  potential 
in  the  subsonic  part  of  the  flow  and  corrections  to  the  shock 
location.  This  idea  will  be  developed  in  more  detail  in  a 
report  where  the  possibilities  of  a  separate  treatment  of 
the  subsonic  and  the  supersonic  regions  with  subsequent  matching 
are  discussed.  The  development  sketched  last  applies  only  to 
that  part  of  the  shock  where  the  flow  behind  the  shock  is 
subsonic. 

The  flow  field  in  the  vicinity  of  the  junction  between 
the  shock  and  the  sonic  line  is  rather  complicated.  The 
author  considers  the  description  developed  by  Nocilla  as 
correct  (Ref.  6)  while  he  believes  that  the  solutions  of 
Germain  (Ref.  7)  which  are  sometimes  quoted  in  this  context 
refer  to  a  different  setting,  namely  the  reflection  of  a 
sonic  line  in  the  form  of  a  shock  of  a  singularity  propagating 
toward  the  sonic  line  (see  the  discussions  in  Ref.  8).  These 
details  are  probably  not  within  the  capabilities  of  present 
flow  computation  methods.  On  the  other  hand,  it  is  questionable 
whether  a  detailed  computation  is  worth  the  effort  from  a 
practical  point  of  view. 

The  determination  of  the  flow  field  described  here  is 
best  carried  out  in  the  form  of  a  Newton  Raphson  iteration 
in  which  one  solves  the  corrections  to  an  existing  approxi¬ 
mation  by  direct  elimination.  In  this  context  the  follow¬ 
ing  observation  may  be  of  interest.  Using  an  extremum  for¬ 
mulation  and  bilinear  elements,  the  author  has  tried  to  solve 
a  purely  subsonic  problem  by  column  relaxation  similar  to 
that  of  Murman  for  finite  difference  equations.  Since  one 
seeks  the  extremum  of  a  certain  functional,  convergence  is 
guaranteed,  just  as  it  is  in  Murman* s  treatment.  However, 
convergence  turned  out  to  be  disappointingly  slow.  This  can 
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For  column  relaxation  one  then  deals  with  a  matrix 

r-8  i 


(1/3) 


1  -8  1 
1  -8  1 


Here  one  has  much  stronger  diagonal  dominance  than  for  finite* 
differences;  that  is,  one  is  closer  to  point  (over)  relaxation. 
This  probably  explains  the  slowness  of  convergence.  The 
situation  is  aggraved  if  one  takes  Ay  >  Ax.  Naturally,  this 
argument  does  not  apply  if  one  uses  direct  elimination 
A  question  of  convergence  still  remains  but  only  because  the 
problem  is  nonlinear. 
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